(N 

o 

O 



o 

I 

H 



(N 
> 

00 

in 

en 

o 

(N 

> 

X 



An efficient iterative method to reduce eccentricity in numerical-relativity simulations of 
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We present a new iterative method to reduce eccentricity in black-hole-binary simulations. Given a good 
first estimate of low-eccentricity starting momenta, we evolve puncture initial data for ~ 4 orbits and construct 
improved initial parameters by comparing the inspiral with post-Newtonian calculations. Our method is the first 
to be applied directly to the gravitational-wave (GW) signal, rather than the orbital motion. The GW signal is in 
general less contaminated by gauge effects, which, in moving-puncture simulations, limit orbital-motion-based 
measurements of the eccentricity to an uncertainty of Ae ~ 0.002, making it difficult to reduce the eccentricity 
below this value. Our new method can reach eccentricities below 10~ 3 in one or two iteration steps; we find that 
this is well below the requirements for GW astronomy in the advanced detector era. Our method can be readily 
adapted to any compact-binary simulation with GW emission, including black-hole-binary simulations that use 
alternative approaches, and neutron-star-binary simulations. We also comment on the differences in eccentricity 
estimates based on the strain h, and the Newman-Penrose scalar ^4. 
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I. INTRODUCTION 

There is a large-scale effort underway to produce models 
of the gravitational-wave (GW) signal from the late inspiral, 
merger and ringdown of binary systems of black holes, cal- 
ibrated against large numbers of numerical simulations 111 . 
These models will be essential to locate and interpret black- 
hole-binary GW signals in the data from second-generation 
laser-interferometric detectors, the first of which, Advanced 
LIGO, may be commissioned as early as 2014 [2-4] . The 
most pressing need is for models of binaries that undergo non- 
eccentric inspiral. It is non-trivial to prescribe initial con- 
ditions that lead to non-eccentric inspiral in numerical sim- 
ulations, and to date there is no systematic procedure to do 
this for simulations performed using the "moving-puncture 
method" [5 6|, which is the most common in the field. 

In fully general relativistic numerical simulations the bi- 
nary's eccentricity cannot be prescribed. The best we can do 
is to use some model to guess initial parameters that may lead 
to low-eccentricity inspiral and, if necessary, adjust those pa- 
rameters until the eccentricity falls below some acceptable tol- 
erance. The problem is further complicated by the difficulty 
of eccentricity measurement. There is no rigorous or unique 
definition of eccentricity for binary systems in general rela- 
tivity. We can define a number of quantities that all reduce 
to the Newtonian limit, and all agree at zero eccentricity [7|, 
but many of these depend on the motion of the black holes, 
which is gauge dependent. (It should be emphasized that all 
coordinate black-hole motion in these simulations is entirely 
due to the gauge variables) It would be preferable to use the 
GW signal, which is far less gauge-dependent. 

Various techniques have been proposed to obtain momenta 
leading to low eccentricity by employing Newtonian or post- 
Newtonian information and short NR simulations Il8l 4l4l . In 
previous work we estimated the initial parameters from solu- 
tions of the post-Newtonian (PN) equations of motion J8]|9). 
For an equal-mass nonspinning binary these resulted in an ec- 



centricity of e ~ 0.0025 (from the NR eccentricity estimator 
that we use in this paper). For larger mass ratios, and for bi- 
naries made up of spinning black holes, the eccentricity was 
larger, even when higher-order PN spin contributions were in- 
cluded [10|. In some cases we further used PN solutions to 
estimate the overall magnitude of the perturbation in the ini- 
tial momenta necessary to correct for the eccentricity iflOl . 
This procedure worked well, but in providing only the magni- 
tude of the momentum adjustment, it was not possible to inde- 
pendently refine both the radial and tangential momenta. The 
method also relied on the gauge-dependent coordinate motion 
of the black holes, which further limited its potential; a second 
iteration of the method was usually not possible, and eccen- 
tricities could not be reduced below e ~ 0.004. 

A powerful iteration method was proposed in Ref. 1(121 . in 
which eccentricities below e ~ 10~ 5 could be achieved in two 
iteration steps. This method was further extended to precess- 
ing binaries in Ref. [14~1 . This method also relies on the co- 
ordinate motion of the black holes. It has been applied to 
simulations that use initial data in particular quasi-equilibrium 
coordinates that are adapted to the motion of the black holes, 
and this means that (a) large non-physical gauge effects in the 
coordinate motion are not apparent, and (b) the phase of the 
black holes' motion can be used from t = 0, which is neces- 
sary for the implementation of the method presented in lfl2l - 
14 1 . These features make it difficult to apply that method to 
moving-puncture simulations. (We discuss this in more detail 
in Sec. |VH| ) In moving-puncture simulations the orbital phase 
cannot be used from t = 0, and we must instead wait until the 
gauge has settled down (at least one orbit into the simulation); 
and even then it contains additional non-physical oscillations 
due to gauge effects. Even in cases where the coordinate mo- 
tion appears to closely reflect the underlying physics, it would 
be preferable to have an eccentricity-reduction method that 
can be applied to the gravitational-wave signal, which is far 
less gauge dependent and is, ultimately, the physically mea- 
surable quantity that we are interested in modeling. 



In this paper we present a robust iterative method that over- 
comes these issues. The idea is as follows. Start with a short 
NR simulation that exhibits eccentricity, and a non-eccentric 
PN/EOB (effective one body, see e.g. lfT5l [T6l) evolution of 
the same system. Adjust the initial momenta in the PN/EOB 
evolution until it exhibits eccentricity oscillations that agree 
with those in the NR waveforms, in both amplitude and phase. 
The inverse adjustment is then applied to the NR initial mo- 
menta, and a new NR simulation performed, and the process 
repeated. The use of the amplitude ami phase of the eccentric- 
ity oscillations makes it possible to independently determine 
the required adjustment in both the tangential and radial mo- 
menta of the black holes. The problem of matching the ampli- 
tude and phase of the eccentricity oscillations can be cast as a 
minimization problem, and its solution semi-automated. 

We describe our approach in more detail in Sec. [TTJ and il- 
lustrate it with simple PN examples (which avoid troublesome 



gauge and noise issues). In Sec. Ill we turn to full NR Simula 



tions. We first describe a method to filter the GW signal, and 
make the crucial observation that the eccentricities measured 
from the phase of the GW strain h and the Newman-Penrose 
scalar 1*4 are not the same; this point is elaborated further 
with a 1PN calculation in Appendix |C 2 1 We then apply our 
method to three NR configurations. 



In Sec. IV we develop a systematic procedure to determine 
the momentum adjustment factors, which makes it possible to 
semi-automate our method. We also make some estimates of 
the computational overhead of applying our method in large 
parameter studies. 

The method could in principle use the orbital phase rather 
than the GW phase. We show in Sec.[V] however, that the or- 
bital frequency contains additional oscillations due to gauge 
effects, which make it difficult to use it for eccentricity reduc- 
tion below e ~ 0.002. 

Our method can reach eccentricities below 10~ 3 in one 
or two steps. In Sec. VI we demonstrate that, perhaps sur- 



prisingly, NR simulations with eccentricities even as high as 
e ~ 0.01 are unlikely to introduce noticeable errors into GW 
searches or parameter estimation in the advanced-detector era, 
and that a target eccentricity of e ~ 10~ 3 reduces phase oscil- 
lations to well below our most stringent current requirements 
on NR phase accuracy. 



II. DESCRIPTION OF THE ECCENTRICITY REDUCTION 
METHOD 

A. Sketch of the eccentricity algorithm algorithm 

In our numerical simulations we start with two black holes 
with masses m\ and mi and spins Si and S2, separated by a 
coordinate distance D. In this work the spins are both parallel 
or anti-parallel to the orbital angular momentum of the binary, 
so there is no precession and the orbital plane is fixed. Given 
such a configuration, our goal is to estimate values of the ra- 
dial and tangential momenta (p r ,Pt) that lead to inspiral of 
(in principle) arbitrarily low eccentricity, i.e. quasi-circular 
(QC) inspiral. In this work we will formulate our method for 



moving -puncture evolutions of Bowen- York-puncture initial 
data, but it can be generalized to other approaches. 

Our method is based on having a sufficiently accurate ap- 
proximate model of the frequency evolution of the GW sig- 
nal as a function of the initial momenta, GM{pt-,Pr), for the 
same initial configuration as used in a numerical simula- 
tion. We choose initial momenta (p%pf) for a first numeri- 
cal simulation, such that the eccentricity in our model is zero, 
e M(p®,P®) = 0. Since the model is approximate, the eccen- 
tricity in the waveform that results from a numerical simula- 
tion using these parameters, e^, R , will be in general non-zero. 
However, we assume that the model, although not precisely 
faithful to a full numerical simulation, does capture much of 
the dependence of the GW signal (and its eccentricity) on the 
initial parameters. 

We then try to remove the eccentricity in our simulations by 
adjusting the initial parameters by the same amount as is re- 
quired to produce the same eccentricity in the model solution. 
This basic idea was already presented in previous work [8 10] 
and we will justify that this is a valid assumption in Sec.|HB| 



In those applications we adjusted the tangential and radial 
momenta by the same factor, i.e., we found a factor X such 
that eM(Xp%Xp®) = e^, R , and then updated the parameters by 

iPriPt) = (Pr ') Pt) I "^ ■ This procedure does not allow for the 
separate identification of p r and p t , placing a lower limit on 
the eccentricity that can be obtained. In addition, we mea- 
sured the eccentricity using the puncture motion of the two 
black holes. This motion is gauge dependent (and indeed the 
motion is entirely due to the gauge choice), and our experi- 
ments show (see Sec.[V]) that this means that we cannot reduce 
the eccentricity to lower than about e ~ 0.002. 

Eccentricity is only uniquely defined for conservative New- 
tonian dynamics. Based on an expansion of analytic solu- 
tions to the Kepler problem for small eccentricities one can 
define eccentricity estimators (see appendix |B l| or [7 1 for def- 
initions). Due to the lack of a unique eccentricity for BH- 
binary evolution it is important to understand how different 



eccentricity estimators are related. In appendix C 1 we make 
such a comparison for a 1PN binary and the GW signal ob- 
tained from the quadrupole formula. To estimate the eccen- 
tricity from NR data, we employ the eccentricity estimators 
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where 0GW,fit(O an d ^fit(0 are approximations to the non- 
eccentric phase and frequency, obtained via a fit over several 
orbital periods. Q. is the orbital frequency and 0qw is me GW 
phase obtained from the Newman-Penrose scalar ^¥4 at some 
fixed extraction radius r ex - The eccentricity is estimated from 
the amplitude of the eccentricity oscillations. 

To determine separately the radial and tangential momenta, 



we consider the eccentricity function Eqn. (2.2 1, which can 
be written with reference to the GW frequencies from the NR 



simulation and the approximate model, 



cnr(0 



2om(0 



(2.3) 



In practice the frequency from the NR waveform, Onr, will be 
very noisy, and it is more convenient to consider the residual 



&(t) = a)m(t)-a)M(t)- 



(2.4) 



If we perturb the momenta by the factors (X,-,Xt), then we can 
also calculate a residual between the perturbed model and the 
non-eccentric model, 

m^{t):=M U {XrMt) = (OM{'Kp rMht)-Ohi{t). (2.5) 

Our modified method then consists of choosing (X r ,X t ) such 
that 



'M 



(0 «#(*). 



(2.6) 



with agreement both in the amplitude and phase of the residu- 
als. Having determined these factors (A^, A r °) for the first NR 
residual &°(t), we produce updated parameters for the next 
numerical simulation, 



pI = p°M 
v\ = rfA°- 



(2.7) 
(2.8) 



We then perform a second numerical simulation using 
(pI >Pt)' If our technique works, then the waveform produced 
in this simulation will contain less eccentricity, e^ R < e^ R . 
The entire process is then repeated, and successive updates 
are made, 



„/+! 



PJ +1 



= pWK, 
= p\IK, 



(2.9) 
(2.10) 



until the eccentricity has fallen below some desired threshold. 

This is our eccentricity reduction procedure. What remains 
is to specify the model of the GW phase and frequency evo- 
lution, procedures filter the NR GW signal for analysis, and 
a method to locate the optimal (X r ,X,) parameters. We now 
focus on these issues. 

We expect that the efficiency of any iterative procedure will 
depend on the fidelity of the model (Bm to the results of a nu- 
merical simulation. In our procedure, we use as our model 
solutions of the EOB equations of motion for nonspinning 
point particles, augmented by the highest known PN spin ef- 
fects. Details are given in Appendix [A] We produce the ini- 
tial guess for (p%Pf) using the same procedure as in our past 
work [8 10 1 : we solve the PN/EOB equations of motion with 
a large initial separation (typically ~ 40M), using initial mo- 
menta from PN circular-orbit expressions (where p r = 0), so 
that the eccentricity has essentially reduced to zero through 
radiation reaction by the time the solution reaches D ~ 10M, 
where we are interested in starting a full numerical simulation, 
and at this point we read off the parameters (p®,pf)- 

Since we base our model for the GW signal on the EOB or- 
bital frequency £2eob we ta k e mt0 account the retardation of 



the GW and the relation between the orbital and GW frequen- 
cies. We define the EOB model 0)m and the NR frequency 
0)nr for the GW frequency at a finite extraction radius r ex as 

to NR (f):=co G w(f + '-ex) (2.11) 

®M(PnPut) :=2£1 EO b(/V,/V>0- ( 2 - 12 ) 

It is also in principle possible to apply the same method 
to the orbital frequency of the puncture motion, or to the 
separation of the two punctures. We find that this works 
adequately well if we place only moderate requirements on 
the final eccentricity. In these cases we use Q.M(p r ,Pt'j) — 
&BOB(Pr,Pt;t), or r M!0rb (p r ,p t ;t) = r EOB , orb (p r ,p t ;t), re- 
spectively, instead of (Oy[{pnPt)- An orbital NR frequency 
residual can then be defined as 



= 2(fl(f)-QM(0). 



(2.13) 



^orb(0 

to replace S% in Eqn. ( 

PN/EOB example II C where we use only orbital quantities, 

and in section III D as a comparison to Si,. In full NR exam- 



2.6b. This residual will be used in the 



pies, however, we will see in Sec. Ill D that this method does 
not allow us to achieve eccentricities below about e ~ 0.002. 



B. Requirements on the model 

We return to the main idea that the eccentricity can be re- 
duced by finding scale factors X that fulfill 



e M (X r p r ,X t p'})=e% R >0 



(2.14) 



and then using [p® / X r , p® / Xt) as new initial parameters. We 
would like to make two points related to this assumption. 

Given two approximants A and B, which could be EOB, PN 
or NR, let one of the two, say A, play the role of the model 
and the other the role of 'NR', i.e., we want to find X such 
that eA(X r p®,X,p®) = e\. Define QC parameters pa,Pb f° r 
each model that satisfy 



eA(pA) =e B {pB) =0. 



(2.15) 



Assuming that there is only a single eccentricity minimum and 
Pa j^ Pb, it then follows that 6a(pb) > and e B {pA) > 0. The 
values eA{ps) an d e B{PA) can be interpreted as a distance be- 
tween the two QC parameter sets pa , ps for the two approxi- 
mants A and B. If these approximants model the binary evo- 
lution up to a similar order of accuracy, then it is natural to 
assume the symmetry 



eA(pB) =e B (p A ). 



(2.16) 



PN and EOB evolutions are symmetric in this sense up 
to numerical accuracy in determining the eccentricity. For 
equal-mass non-spinning inspiral PN/EOB QC data at a sep- 
aration of D = 12M we find that e PN (p EOB ) = 0.003192, 
<?eob(ppn) = 0.003179, an asymmetry of merely 0.4%. This 
distance will vary depending on the location in the BBH pa- 
rameter space, i.e., on the mass ratio and spins of the black 
holes, and on the initial separation. 



The symmetry is weaker between NR and PN/EOB. This 
is to be expected, for two main reasons. The first is physi- 
cal: PN and EOB evolutions differ from each other only in 
higher-order PN contributions, while the NR waveforms cap- 
ture the full general relativistic physics. The second reason 
is related to gauge: the PN/EOB parameters formally map 
to the Bowen- York-puncture parameters only up to 2PN or- 
der ifTTI . and that is only true for the initial data. From the 
construction of the initial data, through the gauge changes that 
the wormholes undergo at early times in a moving -puncture 
simulation [ 1 8—20 1 , up to the point where the gauge has set- 
tled down after approximately one orbit, there do not exist 
any quantitative predictions of the relationship between the 
momenta in a PN/EOB calculation and the physical momenta 
during the inspiral in the moving-puncture simulation. All we 
have are observations that suggest that there is a close rela- 
tionship between the PN/EOB and NR momenta (see Sec. V 
D of Ref. [9]). However, the condition ca{pb) — zb(pa) is not 
strictly required by our method. 

Rather, the crucial assumption is that the behavior of the 
model and its eccentricity Cm(p) in the vicinity of the model 
QC parameters pm is close to that of eNR(p) near its QC pa- 
rameters /?nr- That is, the gradient V p e^ should be close to 
VpeNR in a region around the respective QC solutions extend- 
ing to the highest required eccentricity, say 10~ 2 . We have 
checked that this is indeed the case by comparing -^- for NR 
and EOB data and the explicit Newtonian eccentricity formula 
in the sensitivity analysis given in section |VIA| details are 



given in Sec. VI 



In addition, the model must be sufficiently faithful to the 
real physics to produce reasonable starting momenta for our 
procedure. This rules out Newtonian or 1PN models for our 
method. Obviously, the better the starting momenta, the less 
work is required in reducing eccentricity. Therefore it makes 
sense to use the highest order PN/EOB equations of motion 
available. 



C. An example with two post-Newtonian approximants 

We will illustrate the procedure with a simplified example, 
where a PN solution plays the role of the NR simulation. In 
this way issues of numerical noise and gauge effects are re- 
moved, and we can focus only on the eccentricity reduction al- 
gorithm. The model remains the EOB solution described pre- 
viously. Note that for this illustration, we could equally well 
swap the roles of the EOB and PN solutions due to the symme- 



try discussed in Sec. II B The configuration is an equal-mass 
non-spinning binary with an initial separation of D = 12M. 

For the "NR" simulation using (p®,p®), the eccentricity is 
e « R rv 0.003. The "NR" frequency a>o and its residual 38° are 
shown by the black solid line in Figs. [T]and[2] For the first 
eccentricity reduction step we choose X r = 1 (i.e., we do not 
alter the radial momenta). We make a guess for the pertur- 
bation factor Xt and calculate M^ using the perturbed initial 
parameters. The perturbation factor X, is then adjusted until 



good agreement is achieved. In Section IV A we will describe 
an algorithm to automatize this adjustment; here we will find 



a good agreement between M^ and 38' "by eye". Figs. Ill and 
blillustrate the effect on Q. and M^ of choosing A, = 1.0015 
(green dashed lines). Having obtained acceptable scale fac- 
tors X* = (1, 1.0015), the new initial parameters (p},pj) are 
calculated according to Eqns. j2.1\ and ( |2.8| l. This first step 
results in a reduction of the eccentricity by a factor of 40, and 
the corresponding "NR" residual & 1 is shown by the thick red 
line in Fig. [2] 
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FIG. 1 : GW frequencies in the eccentricity reduction example. The 
grey line indicates the frequency of the non-eccentric reference EOB 
solution, 0>m. The frequency of the surrogate eccentric "NR" simu- 
lation, Wo, is in black. The results of perturbing the initial momenta 
of the reference EOB solution are also shown for X, = 1.0015 (green 
dashed: cofa)- 
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FIG. 2: Same as Fig. H\ but now showing the frequency residuals 
l |2.13fr . The additional red curve shows the "NR" evolution with the 
improved initial parameters. In the first iteration the eccentricity has 
been reduced by a factor of 40, e\ ~ 8 x 10 . 



In the second step of the example, we have to deal with the 
problem of 'dephasing' between the "NR" and model residu- 
als. In Fig. Elwe can already see that the EOB residual (the 
green dashed line) is slightly out of phase with ^°. This de- 
phasing becomes more pronounced as the eccentricity is re- 
duced (see the residual 3& ), and has to be removed if we are 
to continue with the procedure. 



A generic (frequency) residual 

^(f)=Acos(£V + *) + 0(f), 



(2.17) 



is composed of sinusoidal oscillations of frequency £l r and 
amplitude A (which is directly related to eccentricity via equa- 
tions ( |2.3| i) and ( |2.4[ ), plus a non-oscillatory contribution, 
which stems from the different phase evolutions between the 
simulation or perturbed model and the model waveform. It is 
this non-oscillatory contribution that gives rise to the dephas- 
ing in the residuals. 

To resolve this problem, we note that our eccentricity re- 
duction method aims to capture only the eccentricity-related 
oscillations in the residual, and the non-oscillatory part con- 
tains no useful information. To remove the non-oscillatory 
part we perform a fit in time to each residual and subtract it to 
obtain the 'residual modulo dephasing' 



M{t):=m{t)-M & {t). 



(2.18) 



As the fitting model we choose either a polynomial of order n 



^ fit (o=£<v, 



(2.19) 



;=o 



or the rational model given in equation ( |B27| i. A polynomial 
fit is more robust, but its order (usually 4 or 5) needs to be 
adjusted according to the length of the time interval (the 'fit- 
ting window') used for the fit, to avoid picking up parts of the 
eccentricity oscillations. In addition it is advantageous to dis- 
card data in the resulting function 2%(t ) near the boundaries of 
the fitting window to reduce artifacts. 

We now return to the second reduction step in our exam- 
ple. In the top panel of Fig. [3] the solid red line indicates 
the same S% x as in Fig. When we remove the dephasing, 
we recover the solid red line in the lower panel. The eccen- 
tricity oscillations are now clearly visible, and we can again 
search for appropriate perturbations (A r ,A f ) to the reference 
EOB solution to model this residual. We find by trial and 
error that the optimal perturbation parameters are given by 
A* = (A;,A r *) = (1.013,1.000015). Fig.[i]shows the results 
of applying A * to the initial parameters of the reference EOB 
solution. The perturbation ,^Jj matches very well with the PN 
residual. (This perturbation is also shown in the top panel of 
Fig. [3] where we see that the matching cannot be performed 
without first removing the dephasing in the residual £$ l .) Im- 
proved initial parameters are then obtained by adjusting the 
momenta: p r —¥ p r /l.Q13, p t — > p f /l. 000015. These mo- 
menta result in a further reduction of an order of magnitude 
of the eccentricity, which leaves us with ei ~ 8 x 10~ 6 . The 
momenta, eccentricities and scale factors for each iteration are 
given in Tab. II] 
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FIG. 3: Frequency residuals in the second step of the eccentricity 
reduction example. Top panel: Raw frequency residuals 3% (as in 
Fig. 12b, and 3£ 2 , which is the result of this eccentricity reduction 
step. Also shown are the perturbed-EOB residuals, which cannot be 
easily compared to 8$^ . Lower panel: Now the dephasing has been 
removed from 2$} , making it possible to determine the appropriate 
perturbation factors. See text for more details. 



Iteration 



A 



en 



0.000541 0.0851657 0.003 

1 0.000541 0.0850382 8 x 10" 

2 0.000534 0.0850369 8 x 10" 



A/- 



X, 



1 1.0015 
1.013 1.000015 



TABLE I: Initial momenta, eccentricity estimates eo. and results for 
the equal-mass non-spinning PN/EOB example eccentricity reduc- 
tion case discussed in the text. 



produce a clean GW signal, which is the key ingredient in our 
procedure. The procedure itself is then applied to three non- 
precessing black-hole-binary configurations in Sec. |III D| 



A. NR Setup 



III. APPLICATION TO NR SIMULATIONS 

We now apply our eccentricity reduction procedure to full 
NR simulations. In Sec. IIII Al we summarize our numerical 



methods, and in Sec. IIIB we summarize our procedure to 



Our numerical setup is similar to that used in Ref. [ 101, 
but for completeness we repeat the details here. We per- 
formed numerical simulations with the BAM code ll2Tll22l . 
The code starts with black-hole-binary puncture initial data 
generated using a pseudo-spectral elliptic solver [25 1, 



and evolves them with the ^-variant of the moving -puncture 
E version of the BSSN EU |27) formulation of the 3+1 
Einstein evolution equations. Spatial finite-difference deriva- 
tives are sixth-order accurate in the bulk [22 1, Kreiss-Oliger 
dissipation terms converge at fifth order, and a fourth-order 
Runge-Kutta algorithm is used for time evolution. The grav- 
itational waves emitted by the binary are calculated from the 
Newman-Penrose scalar V F4, and the details of our implemen- 
tation of this procedure are given in [21 ]. 

In each simulation, the black-hole punctures are initially a 
coordinate distance D apart, and are placed on the y-axis at 
yi = —qD/{\ +q) and y>2 — D/(l + q), where q — M2/M1 
is the ratio of the black-hole masses in the binary, and we 
always choose M\ < Mi- The masses M; are estimated from 
the Arnowitt-Deser-Misner (ADM) mass at each puncture, ac- 
cording to the method described in [23]. (This measure be- 
comes inaccurate for high spins l28l . but this does not pre- 
clude the application of the eccentricity reduction procedure 
presented in this work.) The Bowen-York punctures are given 
momenta p x = ^fp, tangential to their separation vector, and 
p y = ±p r towards each other. The spin parameter of a BH is 
defined as Xi = St/Mf. 

All simulations used the "X variant" of the moving- 
puncture method and six mesh-refinement buffer points. The 
base configuration consists of l\ nested mesh-refinement 
boxes with a base value of N 3 points, which surround each 
black hole, and h nested boxes with (2A^) 3 points, which sur- 
round the entire binary system. The value of l\ differs for 
each black hole; in equal-mass simulations l\ is the same for 
each black hole, but when the mass ratio is 1:2, the smaller 
black hole is given one extra refinement level, so that both 
black holes are equally well resolved. In addition, we use 
(4JV) 3 points on the level where wave extraction is performed, 
to allow accurate wave extraction to larger radii. The levels 
immediately above and below this are given an intermediate 
number of points (typically (3N) i ), so that no two levels are of 
the same size. The choices of N, l\, h and the resolutions are 
given in Tab. Ill] The resolution around the puncture is denoted 
byMi/h min . 

Far from the sources, the meaningful length scale is the total 
mass of the binary, M = M\+ M%, and so the resolution on the 
coarsest level is given by h max /M. We also give the resolution 
on the wave extraction level(s), h ex /M. 



In section III D we will present results for eccentricity re- 



duction for the configurations listed in table ITTl with the ex- 
ception of the equal-mass non-spinning configuration. The 
latter was used to study the dependence of gauge oscillations 
in the orbital quantities on the parameter tj that appears in the 
T-driver shift condition. These are discussed in Sec. [V] In 
all other equal-mass cases, rj = 2/M, and for unequal-mass 
configurations we used a spatially varying 77 [29-31 1 with the 
functional form 



ri(x) = T] A - 



rip - TlA 



(as?)" 



(3.1) 



fixed the width w = 2.67M and power 8 = 2, so that the mod- 
ification falls off like 1/r 4 . 

Since the focus of this paper is on eccentricity reduction, 
inspiral runs are sufficient. Where available, we also give the 
location of the amplitude maximum in time, t pea k, of the GW 
signal and the number of GW cycles, JVgw for each simula- 
tion. A full convergence series has been performed for the 
q = 2,%\ — 0,^2 = 0.25 configuration only, which is used in 



Sec. VI C to compare phase errors and mismatches due to ec- 
centricity with errors due to numerical resolution. 



B. Filtering the numerical GW signal 

We use the gravitational wave frequency extracted at 
some finite extraction radius as a more gauge-invariant in- 
put quantity. We extract the (£ = 2,m = 2) mode of the 
Newman-Penrose scalar ¥4 and define the wave phase 0gw 
as r*¥f(t) = A(f)e'* cw(f) . The GW frequency 0)gw is then 
obtained as the time derivative of the GW phase. The numer- 
ical noise in 1*4 prevents this definition from being directly 
useful for our setup. This is clear in Fig. [4] where the raw 
GW frequency from ^4, from a simulation of an equal-mass 
nonspinning binary is shown in grey. 

The following fitting and filtering technique has allowed us 
to construct a cleaned residual eccentricity oscillation that can 
be used with our method. We extract the eccentricity oscilla- 
tions from the phase 0gw by first removing the overall behav- 
ior of the phase via a polynomial fit of order n, 

</>Gw.fi,(0 = EC,f" (3-2) 

1=0 

(usually n — 4, 5). Then we smooth the residual phase 

0GW,res(O = 0Gw(O ~ </>GW,fitO) (3.3) 

with a low-pass filter jF/,. We have used either a 'brick-wall' 
filter, where higher frequency modes are simply zeroed out in 
Fourier space, or a wavelet filter using Mathematica's discrete 
wavelet transform with an 8th-order Symlet wavelet (32). We 
then perform a nonlinear fit to a sinusoid, 



^i[0GW, res ](O = Bcos{Q.,-t + 0o) 



(3.4) 



where we have chosen the asymptotic value T]a = 2/M, we 
have set rj p = 3/M near the location x p of the small BH and 



Finally, we take an analytic time derivative of this fit, and 
reassemble a cleaned frequency quantity by adding the time 
derivative of 0GW,fit(O to obtain a cleaned GW frequency, 

fitow.deaned = 0Gw,fit(O - BQ., sin(£2 r f + 0q). (3.5) 

The process is illustrated in Fig.|5]for the NR GW phase from 
an equal-mass non-spinning simulation. Compared to orbital 
quantities, the GW signal is usable only after the junk radi- 
ation has passed, at around 500M. This necessitates slightly 
longer simulations, so that 2-3 periods worth of residual oscil- 
lations are available. The resulting cleaned GW frequency is 
shown along with the raw frequency in Fig. HI It is worth not- 
ing that in some cases we can use the cleaned GW frequency 
at times earlier than the region used for the nonlinear fit, but 
in general this can lead to errors in matching to the phase of 
the oscillations. 



TABLE II: Summary of grid setup for numerical simulations. The grid parameters follow the notation introduced in 1211 : see text. Mi//i m ; n 
denotes the resolution on the finest level with respect to the smallest black hole, while h max /M is the resolution on the coarsest level with 
respect to the total mass, M = M\ +Mi- The outer boundary of the computational domain is at Xj mRX /M, where x, = {x, y, z}- In general l\ 
indicates the number of moving refinement levels around each puncture, and li the number of large refinement levels that encompass both 
punctures. In the q = 2 configurations we use three refinement levels around the puncture of the large black hole, and four around the other. 
hex/M is the resolution on the (main) wave extraction level. The simulations were started at an initial coordinate separation is D and include 
A'gw GW cycles before reaching the amplitude maximum in the GW signal f pea k- The starting momenta are given in sectionllll Dl 
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FIG. 4: The GW frequency from an equal-mass non-spinning simu- 
lation. The black line shows the unfiltered frequency obtained from 
differentiating the GW phase. The red dashed line shows The filtered 
GW frequency Eqn. {331. 
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FIG. 5: The GW phase from the same simulation as in Fig.H] with 
the fit l |3.2fr removed. The blue points show the residual obtained 
after subtracting the fit from the raw NR data. The thick magenta 
line shows the filtered residual. The red dashed line shows the cosine 

fiti^GW.i-eslM- 



Measuring the eccentricity from the GW signal: 
strain h vs ^ 



As mentioned earlier, in our NR simulations we extract the 
GW signal from the Newman-Penrose scalar ^4. We then 
estimate the eccentricity from the phase of the (£ = 2,m = 2) 



mode of ¥4 using Eqn. ( 2. 1 



). 



The measurable quantity in GW experiments is not ^4, but 
the strain h, which is related to ¥4 by two time derivatives. If 
we integrated ^4 twice with respect to time, we could instead 
calculate the eccentricity from the phase of h. 

We might naively expect that the "eccentricity from the 
GW phase" should be the same, irrespective of how the GW 
signal is expressed. But this is not the case. Consider the 
complex GW strain as h — A{t)e'^^> . The eccentricity is re- 
lated to the amplitude of oscillations in <p(t) with respect to 
the phase of a fiducial non-eccentric (quasi-circular) binary, 
(?) = (j>oc{t) +4e*M (f ). If we now consider ¥4 = h and ex- 
press this also as A'{t)e'^ ^'\ then we see that after the applica- 
tion of two time derivatives we will not have 0(f) = '(f). In 
practice the difference between the two phases has been found 
to be very small, but the oscillations in the phase will increase 
in magnitude with each time derivative, and the difference be- 
tween e^m and e^wj will be significant. 

In Appendix |C 2| we consider the GW signal from an eccen- 
tric binary at the first-post-Newtonian (1PN) order, and find 
that the two eccentricities are related by 






1 
4 



3ke 2 



(3.6) 



where k is the periastron advance of the binary, and e 2 = 1 /c 2 
indicates the lPN-order term. This relationship must be borne 
in mind when comparing GW-phase eccentricities calculated 
from the strain and ^4. 

To our knowledge this subtlety of GW-based eccentricity 
measures has not been noted in the literature to date. 



D. Numerical relativity examples 

We are now ready to apply the eccentricity reduction pro- 
cedure to NR simulations. We first present an unequal-mass 
aligned-spin configuration with physical parameters q = 2, 
%\ = 0, and Xi = 0.25. This is the example we will consider 
in the most detail, and in the remainder of the paper we will 
refer to this case as our "reference example". 

As stated in the previous section, we need to skip the initial 
burst of junk radiation and evolve for about 3-4 orbits to get 
accurate estimates of the eccentricity based on the GW signal. 

The PN/EOB QC initial parameters give rise to an initial 
eccentricity (measured from the GW phase) of eo ~ 0.006. 
The top panel of Fig.|6]plots the NR frequency residual cal- 
culated from both the GW phase and the orbital phase. For 
this level of eccentricity, we see that the two residuals agree 
well, and either could be used in the eccentricity reduction 
method. We find that a good match with residuals calcu- 
lated from perturbing the reference EOB solution are obtained 
with k — (1, 1.0028). This match results in the adjustment 
Pt —> Pf/1.0028. The best-match residual ^t^ is indicated by 
a green dashed line. 

The parameters that follow from the first iteration step lead 
to an eccentricity of e\ k, 0.003. The GW- and orbital-phase 
residuals from the subsequent simulation are shown in the 
middle panel of Fig. [6] Varying X to find a good match leads 
to the second adjustment of the initial momenta, /?,-—>■ p r / 1 . 15 
and p t — > p t /0.999. We see in this case that the orbital-phase 
residual now includes some higher harmonics, and its ampli- 
tude is also different to that of the GW-phase residual; it would 
now be difficult to obtain a reliable guess of the perturbation 
parameters based on the orbital-phase residual alone. 

The result of the second iteration step is shown in the lower 
panel of Fig. p] The eccentricity is now e^ ~ 3 x 10~ 4 . In 
principle, we could take another step, and indeed the best- 
match perturbation residual has been calculated. However, 
the eccentricity is already very small, and as we discuss in 
Sec. VI B for eccentricities on that order the uncertainty in 



estimating eccentricity and in calculating a cleaned GW phase 
is very high, around 50%. 

The results for this example case are summarized in 
Tab. 



IUa 



The progress of the eccentricity reduction method 
for two additional aligned-spin configurations is given in 
Tabs. |IIIb| and |HIc| For all cases we provide eccentricity es- 
timates from e^.Gw along with the radial and tangential mo- 
menta and their scaling factors for each reduction step. To 
provide another example of the presence of higher harmon- 
ics in orbital quantities, we also give eccentricity estimates 

and plot the or- 



Illa 



from en(0 for the example case in Tab. 
bital frequency estimator for the iteration steps in Fig. 17] This 
plot again indicates the unsuitability of the orbital phase for 
this method. We can see from both Fig. [JJand Tab. |IIIa| that 
the orbital-motion-based eccentricity estimator does not fall 
below 0.0013, which, if it were correct, would suggest a GW- 
phase eccentricity of around 0.002. This is the basis of our 
claim that the orbital motion cannot be used to reduce to the 
eccentricity to below e ~ 0.002. We will study the behavior 
of the orbital motion further in Sec. [V] 
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FIG. 6: Eccentricity reduction steps for q = 2,X\ = 0,^2 = 0.25. 
At each step the NR residual ,^f" is calculated from the filtered GW 
signal Ogw ( re d, thick line), and for reference we also show !%' olb 
calculated from the orbital frequency SI (thin blue line). The best- 
match perturbation parameters X* lead to the residuals 31%, (green, 
dashed line), and are given in Tab. |IIIa| which also provides the up- 
dated momenta. The eccentricities are 0.006, 0.003, and 3 x 10~ 4 . 



We give two significant digits in 1 — X . This choice is sen- 



sible in light of the discussion of errors in Sec. VI A Lastly, 
the uncertainties in the eccentricity values are about 25% for 
e0,GW and, if we do not incorporate the contribution of gauge 
harmonics in the orbital frequency, about 5 — 10% for en. 
These uncertainties are discussed further in Sec. I VI Bl 

In summary, the data presented in this section demonstrate 




400 



600 800 

t[M] 



1000 



FIG. 7: Eccentricity estimator based on the orbital motion, ea for 
the example configuration. The oscillation period is about 330M and 
the dominant feature in the original and first-iteration data, shown in 
grey orange-dashed lines, respectively. At the second iteration step 
higher harmonics dominate and result in an eccentricity (red, thick, 
small-dashed) that is considerably higher than the result of e^ gw- 
Refer to table llllal for details. 



that starting from EOB QC initial parameters we can in gen- 
eral reach eccentricities well below 10~ 3 for aligned-spin NR 
configurations in two iteration steps. 



Iteration 



/>!■ 



Pt 



e 0,Gw e n 



0.000758 0.11710 0.006 0.0045 

1 0.000758 0.11677 0.003 0.0029 

2 0.000660 0.11689 0.0003 0.0013 



A r Af 



1 1.0028 
1.15 0.999 



(a) Results for the example NR eccentricity reduction case discussed 

in the text, q = 2,%\ = (),Xl = 0.25. We give eccentricity estimates 

based on both «o. and 0qw ■ 



Iteration p r p, 


«0,GW 


A r A/ 


0.000677 0.11466 


0.0033 


1.2 0.9985 


1 0.000562 0.11483 


0.002 


0.9 1.0006 


2 0.000623 0.11476 


~ 0.0008 





(b) Results for the configuration 
q = 2,X\ =-0.75,^ 2 = -0.75. 



Iteration 



Pi 



e fGW 



0.000647 0.08764 0.006 

1 0.000539 0.08737 0.001 

2 0.000415 0.08737 0.0003 



A r A{ 



1.2 1.003 

1.3 1 



(c) Results for the configuration 
q= 1,*! =0.5,22 =0.5. 

TABLE III: Results of the eccentricity reduction method for three 
aligned-spin BBH configurations. We give eccentricity estimates, 
the initial momenta p r ,Pt and the obtained scaling factors X r ,Xt for 
EOB QC parameters and two iteration steps in each case. 



IV. IMPLEMENTATION ISSUES 
A. Determination of the scale factors 

In the previous PN and NR examples we determined the 
momenta scale factors (X r ,Xt) by trial and error. In this sec- 
tion we present a systematic procedure that can be automated. 

Recall that our goal is to "match" the NR and model fre- 
quency residuals modulo dephasing 



■ M-Mu 



%M- 



— ^N/T '^\ 



MfiP 



and Eq. dZ6| translates to the requirement that 



2>X 



(4.1) 
(4.2) 



(4.3) 



In achieving this, the crucial observation is that the radial 
and tangential momenta perturbations make contributions to 
the (oscillatory) frequency residual M M that are out of phase. 
Consider a Newtonian binary with perfectly circular orbits. If 
we perturb the tangential momentum at t = 0, then the result- 
ing binary will now follow an elliptic orbit, and the location 
of the bodies at t = will correspond to an extremum of the 
separation, and also an extremum of the instantaneous orbital 
frequency. This point will therefore be an extremum in the 
frequency residual Mm calculated with respect to the circular 
orbit, and we can write the residual as Mm = A t (X t )cos(Q. r t), 
where Q. r is the frequency of the eccentricity oscillations, and 
A t is the amplitude of the oscillations due to the perturbation 



of the tangential momenta. Similarly, a perturbation of the ra- 
dial momentum will lead to an elliptic orbit in which t = cor- 
responds to an extremum of the radial velocity, and therefore 
a zero in the frequency residual. The frequency residual can 
then be written as Myi = A r (X r ) sin(n,f ). A more detailed and 



qualitative exposition of these points is given in Appendix B 2 



For a general perturbation of the momenta, we may then 
write the residual as 



# M = A r (X r )sm(Q. r t)+A t (X,)cos(Q. r t) 
= Acos(Q. r t + A<P). 



(4.4) 
(4.5) 



We see in Appendix B 1 that the amplitudes of the PN/EOB 
residuals depend linearly on the momentum perturbations, 
i.e., A r — a r (X r — 1) and A t — a t (X t — 1). In principle, then, 
we could measure the amplitude and phase of the numerical 
residual M l , and having determined a r and a,, could analyti- 
cally calculate the appropriate scale factors (X ri Xt). 

This method can work well, but in practice we found that 
a more robust procedure consisted of simply performing two 



line searches: first vary the amplitude A in Eq. (4.5 1, but with 
the phase A4> fixed, and then vary the phase with the ampli- 
tude fixed. The searches work as follows. 

Given the NR and model residuals, M'{t) and Mm{^), we 
pick one extremum in M 1 (t ) in the middle of the available time 
window and record its location f, in time and its amplitude 
value Aj — M l {ti). Next, we locate the closest extremum with 
matching sign in Mm and again record its location fjyi and its 
amplitude Am = My\(ty\). We typically store a list of data for 
one or two extrema to the left and right of the fiducial ones, 
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and calculate the average. We then define the deviations d& 
and d$, as 



d A = 1 -A M /Aj, 
d<$> = (t{ — tu)/T, 



(4.6) 
(4.7) 



where bars denote averages over the chosen extrema and T is 
the average time period of the oscillations. It should be clear 
that dA will be zero when the amplitudes are equal, and c/<j> 
will be zero when the residuals St (t) and Myiif) are in phase. 

We first perform a line search for amplitude adjustment. We 
set X r — 1 and choose a symmetric interval in X, around unity 
(with length on the order of 0.005). We then choose X, smaller 
or larger than unity depending on which choice gives better 
agreement of the phase of the residuals (i.e. where \d$>\ is 
smaller). Since dA is a continuous function of A and we have 
a bracketed root, we can then use a robust root-finding method 
that uses relatively few function evaluations, such as Brent's 
method [33], to quickly align the amplitude of the residuals. 

The optimal amplitude from this first search then serves as 
input for the phase adjustment. By construction A<j> is periodic 
in the phase angle 9 = Arg(A r + iA t ) = atan2(A,,A r ) and is 
only continuous for comparison of a given set of extrema in 
both residuals, but not when passing to the next set of extrema 
in one residual due to adjusting the phase. If we are careful to 
select a bracket in a region of continuity of 6, we can also use 
Brent's method for the phase adjustment and quickly find an 
optimal value for A . In practice some care must be taken to 
choose the size of the parameter intervals for the line searches. 
Otherwise the algorithm is automatic. 

An example for the locations covered by the amplitude and 
phase adjustment line searches using Brent's method in the A- 
plane is given in Fig. [8] In this example X, dominates and A, is 
close to unity, which is typical for eccentricities ~5x 10~ 3 . 
This is for the same q — 2,Xi — 0,Xi — 0.25 configuration 



as discussed in Sec. HID There, the adjustments were ob- 
tained manually, while the automatic search was developed 
later, and so the updated parameters are somewhat different to 



those given in Tab. Ilia 



The eccentricity reduction algorithm described here has 
been implemented in Mathematica and will be made available 



at http://gw-models.org The automated method to determine 
the optimal momentum scale factors {A* , X*} requires further 
optimization against a wider range of binary configurations, 
but will also be made available in due course. 



B. Computational cost of eccentricity reduction 

In this section we make a rough estimation of the relative 
cost of eccentricity reduction, within the overall production of 
high-quality numerical-relativity waveforms. 

Our experience suggests that we can measure the eccentric- 
ity with sufficient accuracy from the lowest resolution simu- 
lation of a convergence series. Furthermore, we assume that 
we will carry out two iteration steps, each 3-4 orbits long. As- 
suming that the full simulations will cover approximately 10 
orbits l34l[35l . the length of these two simulations together 




FIG. 8: Line searches in the plane of perturbation scale factors 
(A r ,A/) for the first iteration of the q — 2,%\ = 0,#2 = 0-25 config- 
uration. An amplitude search to minimize d&(X) (blue spheres with 
orange filling) and a phase search to minimize d$(X) along a line 
of constant amplitude (magenta spheres) are carried out. The final 
result of the two searches are the scaling factors A = {1.12, 1.0024} 
(larger red sphere). 



shall be comparable to the merger time for the configuration. 
For simplicity, let us choose the scale factor between the grid 
resolutions used in the convergence series about p ~ 1.15, 
although it usually varies slightly between different resolu- 
tions. The computational cost for each resolution is then pro- 
portional to (l,p,p 2 ) 4 , about (17%, 30%, 53%) of the total 
cost. A conservative value for the relative cost overhead from 
eccentricity reduction is then about 20%. 

For clustered configurations in a parameter study it can be 
argued that instead of a full convergence series a single high 
resolution simulation is sufficient for a subset of these con- 
figurations. In that case the overhead cost for eccentricity re- 
duction is higher, about 35%. On the other hand, for longer 
waveforms or higher accuracy requirements on the final wave- 
form, the overhead cost will be even lower. More importantly, 
in unexplored parts of the parameter space, it is usually neces- 
sary to perform test runs before launching production simula- 
tions. The cost of eccentricity reduction can then be partially 
absorbed into such test runs. This situation will be common, 
and therefore we estimate the overhead cost from eccentricity 
reduction in a large parameter study as between a quarter and 
a third of the total computational cost. 



V. GAUGE DEPENDENCE OF THE ORBITAL MOTION 



In the NR examples in Sec. HI D| we saw that the orbital mo- 
tion cannot be used to make a reliable estimate of the binary's 
eccentricity. This is clear in Figs. [6] and H\ 

In this section we explore in more detail the effects of the 
gauge choice on the puncture motion. The puncture motion is 
generated by the T-driver condition [36], which, in its usual 
formulation, has one free parameter tj. We study the accu- 
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racy of the orbital frequency eccentricity estimator e^, as a 
function of standard choices of tj . 

We first establish the theoretical relationship between the 
eccentricities calculated from orbital and GW quantities. 

The fundamental frequency of the oscillations due to resid- 
ual eccentricity, Q. r , is related to the average orbital fre- 
quency £2<j> via the fractional periastron advance per orbit 
k = A<S>/(2n) by Eqn. (|C9) (see also Q), 



Or 



= l+jfc(ftr). 



(5.1) 



As a consequence, eccentricity estimators for the orbital 
phase and frequency, or GW phase and frequency, are also 
related via this ratio, 






l+k. 



(5.2) 



However, the ratio between eccentricities calculated from the 
orbital frequency and GW phase (i.e., one quantity is from 
the orbital motion and the other from ^4) is different. In ap- 



pendix C 1 we show that to 1PN order and for low frequencies 
(or small k), the ratio is 



k:= 






21 
16 



32 4,' /a 



(5.3) 



where v is the symmetric mass-ratio and e = 1/c. For larger k 
it is more accurate to form the ratio directly from Eqns. \C\A\ 
and( |C24") . 

We now consider a series of equal-mass non-spinning con- 
figurations with spatially constant tj ranging from zero to our 
standard value of TJ = 2/M. Simulations were performed for 
two different sets of initial parameters, one which lead to mod- 
erate eccentricity e ~ 0.004 (Fig. 15), and another for very low 
eccentricity e ~ 0.0005 (Fig.[T(|. 

To reliably assess the residual eccentricity we have com- 
pared eccentricities computed from en and e^.GW for the same 
fitting window t E [400, 1200]M. e<j, gw nas been computed 
for extraction radii R ex = {40,50,60,80,90}M, while taking 
into account the retardation of the wave signal. We have 
computed the orbital eccentricities as the average between the 
minimum and the maximum of the estimators in the windows, 
and have made the required quasi-circular fits in Eqs. ( |2.1[ ) and 
(2.2 1 using a fifth order polynomial for e^ow and a fourth- 
order polynomial for e^. Additionally, e^cwt/) was esti- 



mated using the smoothing procedures discussed in Sec. Ill B 



The average orbital frequency for the NR q = 1 example 
in the fitting window is MQ. s=s 0.027. Taking the fractional 
periastron advance, k, from the Fig. 7 of Ref. [7], a value k ~ 
0.38 seems reasonable. This gives a factor K « 1.25. (For 
comparison, the 1PN definition gives k sa 0.23 which leads to 
1.27). This value of K has been used in Figs.[5Jand 10 



to 



scale ea so that a direct comparison with e^ow is possible. 

Convergence in e^.GW as a function of the extraction radius 
for fixed 77 is spoilt by uncertainties introduced by the level 
of noise in the waveform at these small eccentricities and by 
uncertainties from the sinusoidal fits to the data. 



The resulting eccentricities in Fig. 15] are consistent within 
10% for both estimators and the dependence on tj is very weak 
at this moderate eccentricity. In Fig. [TO] the rescaled orbital 
eccentricities align well with e^GW near f] ~ 0.5. For higher 
values of 77 towards the standard value 2/M the amplitude of 
the second harmonic increases (see also Fig. [TT) . For 77 = 
all quantities are extremely noisy, and en(t) had to be filtered 
and it was not possible to compute e^ gw- 

We find that the absolute error e^ is roughly constant be- 
tween the two cases; it is simply far more noticeable in 
Fig. [TO] where the eccentricity is roughly an order of mag- 
nitude lower than in Fig. |5J 

To see which frequency components contribute to the ec- 
centricity values of en, Fig. [TT] shows the ten lowest Fourier 
amplitudes in en(f) for the simulations shown in Fig 
a function of tj. The amplitude of the harmonic at < 
increases with tj and for tj 
mental frequency. 
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FIG. 9: The eccentricity estimators eQ and e$ gw (f° r extraction radii 
Rex = 40, 50, 60, 80, 90M) as a function of the gauge parameter 77 
for equal-mass non-spinning evolutions with moderate eccentricity 
initial data, cq has been rescaled by a factor K = 1.25 (see text). 
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FIG. 10: Same as Fig. [9] but for an eccentricity almost one order of 
magnitude lower. 
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if A,- = 1, the Newtonian eccentricity formula ( |B17| i gives 
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FIG. 11: The amplitude of the lowest Fourier modes of the eccentric- 
ity residual ea(?) for equal-mass non-spinning evolutions (M = 1) 
with varying gauge parameter r\ . The average orbital frequency is at 
Mco = 0.027. 



One implication of the above analysis is that eccentricity 
estimators computed from the orbital motion become unreli- 
able below a certain eccentricity, e s=a 0.002, and will lead to 
estimates higher than the physical eccentricity present in the 
evolution. Having shown consistency between en and e^ gw, 
we rely in the following on the estimator e0,ow(O- A study 
of Figs. [9] and 10 along with Fig. 1 1 might suggest that the 
orbital-frequency eccentricity estimator could be used if we 
used, for example, 7} = \/M instead of rj = 2/M. However, 
we only know that this value of r\ gave a reasonable estimate 
of the eccentricity after performing a detailed study the 7} de- 
pendence of the estimator for this case. The optimal choice 
of r\ for another configuration could be quite different, and 
may change yet again when a spatially varying 77 is used. The 
GW-frequency-based eccentricity estimator is clearly far less 
dependent on the choice of gauge (as we would expect), and 
we choose to use that in all subsequent work. 



VI. ERRORS 

In the following sections we discuss sources of error in the 
determination of the scale factors, the measurement of the ec- 
centricity, and in the final waveforms. 



A. Sensitivity of eccentricity on scale factors 



dkt 



■■2kt 



(6.1) 



Let pi denote the QC value of the momentum p,, so that 
e(p°,p°) = 0. In NR cases it is easier to estimate Ae/Ap t « 
16. Using the chain rule and the relation p, = A./?° we have 



de 



de dp, 

dp, d X, 



dp.' 



(6.2) 



and find 4§- « 1.9, which is very close to the Newtonian sensi- 
tivity. (Note that this agreement is the basis of our of statement 
that V p eM s=s V p eNR in Sec. II B ) Due to the lack of radiation 
reaction, the Newtonian model fails to give a useful estimate 
of the sensitivity due to perturbations of the radial momen- 
tum. In order to get a theoretical estimate of this value, we in- 
stead performed a sample of EOB/PN perturbations and found 
4r sa 2 and 4r ~ 0.005 for our reference example. 

We can now estimate a lower bound on the achievable ec- 
centricity with our method, using the EOB/PN sensitivities 
and the fact that the perturbation factors are usually very close 
to unity (typically, 1 — X r ~ 0.1 and 1 — kt ~ 0.001). Assum- 
ing an error of x% in X t , the error in the eccentricity will be 
Ae = §£-AXt ~ Jc/50. For Ae < 0.001, we must therefore have 
x < 0.05%, i.e., X, has to be accurate to 5 x 10~ 4 , and by a sim- 
ilar argument X r has to be accurate to within 15%. Clearly, ec- 
centricity is a lot less sensitive to changes in X r than to Xt, but 
on the other hand, the adjustments are larger for X r . While our 
automatic procedure adjusts both momenta in each step, we 
have found that X, is far more important, especially in the first 
iteration step. This is the reason why X r is sometimes equal to 



unity in the example configurations shown in Tabs. Ilia Illb 



and |IIIc| where the adjustments have been checked by eye. 

This raises the obvious question: what is the minimum 
eccentricity we can obtain with p r — 0? If it were suffi- 
ciently low (for example, less than 10~ 3 ), then our eccen- 
tricity reduction procedure could be dramatically simplified 
by locating only the appropriate tangential momentum. This 
can easily be tested on the PN level or estimated using the 
computed sensitivities and a required Ap r = p r . We calcu- 
late jf- w 8 and have Ae s=s j§-Ap r s=s 0.005 for the reference 

(q =2,%\ = 0,^2 = 0.25) configuration. This is lower than 
the eccentricity that often results from the first-guess EOB/PN 
parameters, but is higher than the 10~ 3 threshold that we have 
been aiming for. We also expect this value to be higher for 
configurations with high anti-aligned spins, where the inspiral 
is more rapid, and therefore the appropriate p r is higher. 



We first consider how the eccentricity depends on the un- 
certainties in the scale factors. 

From Newtonian considerations one can see explicitly that 
the variation in the eccentricity depends linearly on the scale 
factors, i.e., 4$- is constant for e ~ 0, where X. denotes either 
of X r or X,. The absolute error in the eccentricity is then pro- 
portional to the error in A.: Ae = ^-AA. °< AX,. For example, 



B. Finite-extraction-radius and fitting errors 

We extract the gravitational wave signal at a number of fi- 
nite radii. For each radius we can compute an approximate 
match and scale factors (A, ,A r ). How large is the related error 
in the scale factors, compared to their required accuracy in or- 
der to attain a given eccentricity? An analysis of iteration 1 of 
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the reference example shows that this error is about 2.5 x 10~ 4 
in Xt and 0.07 in X r . Comparing with the results from the sen- 
sitivity analysis above, this level of error will allow us to attain 
an eccentricity on the order of e ~ 5 x 10~ 4 , although in prac- 



tice we seem to do a little better (see table Ilia I. At this point 
it is difficult to say whether numerical noise, errors due to fit- 
ting artifacts, or errors due to finite extraction radius are the 
dominant source that limits the accuracy of the iteration steps. 
Apart from the errors due to gauge oscillations in ejj, which 
we have pointed out above, eccentricity estimators introduce 
an additional source of error through the polynomial fits in- 
volved. Two relevant parameters are the order of the fitting 
polynomial and the duration of the fitting window. We choose 
the order of the polynomial based on the signal length (2-3 
periods) in order to avoid also picking up the eccentricity os- 
cillations. For our preferred estimator e^cw we also apply a 
low-pass filter to the residual, perform a nonlinear fit to a si- 
nusoid, and take its amplitude as the eccentricity value. From 
tests with a NR signal with an artificially injected eccentric- 
ity (i.e., a sinusoidal modification to the GW phase), we find 
that e^.GW tends to underestimate the eccentricity by 5 — 20%. 
Since real NR data will have less than perfect sinusoidal oscil- 
lations, we take a conservative error estimate of 25% relative 
error in e^.Gw- F° r very low eccentricities c~5x 10~ 4 , the 
amount of noise in the NR phase degrades the accuracy further 
and we estimate the error to be about 50%. 



C. Phase errors and mismatches between eccentric hybrids 

The motivation for this work is waveform modeling for 
gravitational-wave astronomy. Most black-hole-binaries vis- 
ible to the second-generation ground-based detectors Ad- 
vanced LIGO and Virgo will follow non-eccentric orbits. We 
therefore want to simulate non-eccentric binaries. Since we 
cannot construct simulations of binaries with precisely zero 
eccentricity, we are then faced with the question: what level 
of eccentricity in NR waveforms is tolerable! More precisely, 
if NR waveforms of a given eccentricity are used to produce 
waveform models, which are in turn used for GW searches 
and parameter estimation, how many signals will the eccen- 
tricity cause us to lose, and of those signals that we observe, 
how much will our parameter measurement be skewed? (The 
effect on detection of the opposite problem — using non- 
eccentric models when the real signals are from eccentric bi- 
naries — is considered in [37 1.) 

The standard tool with which to address such questions is 
the mismatch. We choose one waveform, hi, as the true signal, 
and another h 2 as the model that will be used in GW searches 
and parameter estimation. The most basic mismatch is defined 
as 

(hi\h 2 ) 



= 1 — max - , 



(6.3) 



where the inner product between the two waveforms is further 
defined as 



The inner product is calculated in terms of the frequency- 
domain waveforms h(f), and is weighted by the power spec- 
tral density S„ (/) of a given detector. The frequency range in 
which the detector is deemed sensitive is [/ m ; n , /max] ■ In calcu- 
lating the mismatch, the inner product is maximized over time 
and phase offsets of the waveforms, so that two waveforms 
produced by exactly the same source (even at different times 
and with different initial phases) would lead to a mismatch 
of zero. In our calculations we use for S„ (f) the Advanced 
LIGO zero-detuned, high-power |[39ll noise curve, and choose 
/min = 20 Hz, and / max = 8 kHz. 

In assessing the usefulness of a given model in a GW 
search, we should also maximize over the physical parameters 
of the model. By not doing that, we produce an upper limit on 
the mismatch between the signal and the model. We do not 
perform the additional maximization because (a) we cannot, 
because our NR waveforms only model discrete configura- 
tions, and (b) we can also use the mismatch ( |6.3[ > to estimate 
the indistinguishability of the waveforms, which characterizes 
the effect of the waveform's error on parameter estimation. If 
h\ is observed at a given signal-to-noise ratio (SNR), then it 
cannot be distinguished from h 2 if |5/z| = \h\ — h 2 \ 2 < 1 ||40l . 
If this is the case, then the parameter estimation uncertain- 
ties will not be effected by any errors in the model hi. The 
indistinguishability can be related to the mismatch by fiTll 
\8h\/p 2 = 2Jt, where p is the SNR. Note that for any h 2 
there will be an SNR sufficiently low that we cannot distin- 
guish it from hi, and an SNR sufficiently high that we can 
distinguish them (unless of course they are in fact identical). 

The primary application of our waveforms will be in pro- 
ducing phenomenological waveform models R214461 . These 
are based on hybrids between PN early-inspiral waveforms 
and the NR late-inspiral-and-merger waveforms. In calculat- 
ing mismatches, we will therefore consider PN-NR hybrids of 
our waveforms. And since we are aiming for a non-eccentric 
model, the PN portion of the waveform will always have zero 
eccentricity. We produce hybrid waveforms using the stan- 
dard technique where we align the PN and NR waveforms 
over a small time interval around the frequency Mco — 0.055 
and then smoothly blend them together. We consider the ref- 
erence example (q = 2,^i = 0,^2 = 0.25) configuration, and 
use TaylorTl as the PN approximant. 

Fig. 12 shows the mismatch results between a "model" 

= 0.01, and a "signal" wave- 
For comparison we 
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hi(f)h* 2 (f) 

Sn(f) 



df 



(6.4) 



waveform with NR eccentricity e - 
form with NR eccentricity e = 3 x 10~ 4 . 
also show the mismatch between simulations with different 
levels of numerical resolution. The precise mismatch values 
in the figure should not be taken too seriously; these are sen- 
sitive to the details of the hybrid construction and the mis- 
match calculation. The main point is that the mismatches are 
remarkably low. For reference: (a) a commonly used crite- 
rion for detection requires mismatches between the signal and 
model of less than 0.03; (b) the mismatches between a series 
of equal-mass nonspinning-binary NR waveforms, each pro- 
duced with a different NR code, were found to be typically 
five times larger than the 10~ 4 that we see here (47), and mis- 
matches for larger mass ratios and spinning black holes tend 
to be higher [46, 48 1; and (c) these waveforms would be in- 
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distinguishable for SNRs below 100 in Advanced LIGO. An 
observation with such a high SNR in Advanced LIGO would 
be truly exceptional. In third-generation detectors, such as the 
Einstein Telescope (ET) [49], and space-based detectors (for 
example |50|), however, far larger SNRs would be likely, and 
these differences will be distinguishable. 
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FIG. 12: We show a mismatch between PN-NR hybrids of ec- 
centricity e = 0.01 and a reference hybrid (e ~ 0.0003) for the 
q = 2,%\ =0,X2 = 0.25 configuration (red curve). The hybrids use 
the same numerical resolution (N = 88 grid) and extraction radius. 
For the sake of comparison, a second mismatch is computed between 
different resolutions N = 88 and N = 96 for eccentricity e ~ 0.0003. 

This plot suggests, then, that the use of NR waveforms with 
eccentricities as high as even e ~ 0.01 will have negligible 
impact on GW detection and parameter estimation over the 
next decade. So long as the raw PN/EOB momenta parameters 
lead to eccentricities below this value, we could conclude that 
we do not need to apply any eccentricity reduction at all. 

The weakness of mismatch plots is that they are effectively 
a measure of how well the phases of two waveforms can be 
aligned, and no more. In particular, they do not tell us how 
the errors in a set of waveforms will affect the calibration of 
coefficients in a phenomenological or EOB model. This ques- 
tion cannot be answered until we attempt to calibrate those 
models. What we do know is that the eccentricity reduction 
process is relatively computationally cheap to perform (see 



Sec. IV B I, but repeating a large family of simulations if they 
turn out to be inadequate is extremely computationally expen- 
sive. For this reason, we take the conservative view that we 
want the phase errors due to eccentricity to be no larger (and 
preferably smaller) than those due to other numerical errors. 

This requirement motivates Fig. [13] in which we compare 
the v ?4 phase differences between inspiral NR waveforms of 
varying eccentricity. All curves are with reference to the low- 
est eccentricity phase (e~3x 10~ 4 ). The phases are aligned 
at Mco = 0.055 using a non-eccentric fit through the data. We 
see that the phase differences have a secular and an oscillatory 
part. Since the amplitude of the phase oscillations is linearly 
related to the eccentricity, we were also able to estimate the 
curve for a simulation with our threshold eccentricity require- 
ment of e = 10~ 3 (the short-gap-dashed magenta line). 

Note that simple alignment at a fixed frequency can lead 



to large secular phase errors that depend sensitively on the 
alignment frequency. The worst-case phase errors due to ec- 
centricity for this set of simulations agree well with a Caltech- 
Cornell estimate in Eq. 61 in iTPSll . As an example, for phases 
aligned at Mco = 0.1 the backwards-in-time phase error for 
e = 0.003 reaches 0.8rad at Mco = 0.055. But it is impor- 
tant to realize that this is a pessimistic estimate. By using 
non-eccentric fits to the phase evolution, we have effectively 
aligned to an underlying "mean quasi-circular frequency", and 
now the secular drift is far less. This also mimics what is 
effectively done in producing hybrid waveforms, where the 
root-mean-square integrated phase disagreement is minimized 
over a time interval that includes at least one GW cycle. 

Based on our previous experience with moving-puncture 
simulations, we expect the accumulated numerical phase error 
through inspiral to be less than 0.01 rad, and this is shown by 
the shaded region in the figure. We see that the oscillations in 
the phase disagreement are greater than 0.01 rad for the sim- 
ulations with e = 0.003,0.006,0.01, but are well within this 
tolerance f or e = 10~ 3 . 

This figure does not include dephasing through merger and 
ringdown. We expect that the eccentricity will have negligi- 
ble effect on the merger/ringdown waveform, but the phase 
oscillations through inspiral that are visible in Fig. [13] may 
cause the eccentric binary to merge slightly early or later 
than its non-eccentric counterpart. We can estimate this de- 
phasing effect as follows. Based on the figure, we see that 
near merger (we choose Mco = 0.1), the non-eccentric and 
e = 10~ 3 binaries may be as much as 0.01 rad out of phase, or 
~ l/500th of a cycle. The orbital period of the last full orbit 
is roughly T ~ 100M, and so this dephasing corresponds to a 
relative time lag of 8t m Q.2M. If we take the GW frequency 
function for a non-eccentric waveform, co e= o(t), and integrate 
(co e= o(t - 8t) - co(t)) from Mco = 0.1 through merger and 
ringdown, we can estimate the additional accumulated de- 
phasing. This effect will be largest when the ringdown fre- 
quency is highest, i.e., for an equal-mass binary with large 
aligned spins. Using the results from a previous simulation of 
a q = 1, Xi ! = 0.85 configuration [ 10], we find an additional ac- 
cumulated dephasing of 0.2 rad. This is well below the lowest 
accumulated numerical phase error that we have recorded to 
date (1.0 rad for the q= 1 , J, = 0.5 configuration, in Tab. Ill of 
Ref. iflOl ) and suggests that our threshold of e = 10~ 3 limits 
the phase uncertainty due to eccentricity to a level lower or (at 
worst) comparable to the numerical phase error. 

There are two points that should be made about this esti- 
mate. We should recall that this is a conservative requirement 
on the eccentricity; the mismatches shown in Fig. 12 suggest 
that a requirement of e < 10~ 3 is much lower than what is re- 
quired for GW astronomy in the advanced detector era. On the 
other hand, if one does wish the eccentricity-induced phase ef- 
fects to be below numerical uncertainty, then the eccentricity 
threshold will be different for simulations at different levels of 
accuracy. If the numerical phase errors were an order of mag- 
nitude lower (as in the simulation discussed in Ref. [51 1) then 
the same argument would demand e < 10~ 4 , which is what 
was used in that case. 
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FIG. 13: Phase differences from ¥4 for q = 2,%\ = 0,%% = 0.25 
simulations of varying eccentricity with N = 88 numerical grids. All 
differences are with respect to the e* .gw = 0.0003 simulation. We 
also estimate the curve for a simulation with eccentricity e^cw = 
0.001 (magenta, short-dashed); see text. A stringent NR phase error 
requirement of ±0.01 rad is indicated by the shaded region. 



VII. COMPARISON BETWEEN ECCENTRICITY 
REDUCTION METHODS 

We now compare our new techniques with the iterative 
eccentricity-reduction method used by the Caltech-Cornell- 
CITA group in Refs. lfl2l4T4ll . We will refer to this as the 
"CCC" method. The CCC-method was able to efficiently re- 
duce eccentricity down to the order of 10~ 5 for generalized- 
harmonic formulations of Einstein's equations [52 1 using 
conformal-thin-sandwich (CTS) excision initial data ll53l . 
Since this method has proved so successful, it is important 
to explain why we have gone to the trouble of devising an 
entirely new method, and why we have not simply adopted 
their method. The key reasons are associated with the gauge 
choice, and we explain this further below. First we provide a 
brief description of the CCC method. 

We discuss here only the latest version of the method lfl~4l . 
In a nutshell, the time-derivatives of the coordinate separa- 
tion and the orbital frequency are fit against model functions 
consisting of a non-oscillatory part that approximates the in- 
spiral, plus a sinusoidal term to capture the residual eccen- 
tricity oscillations. The fit results are then used to construct 
updating formulas for the initial velocity and the initial orbital 
frequency at a starting separation and partly make use of in- 
formation from the Newtonian Hamiltonian. 

In principle we could use the CCC method with the 
moving -puncture orbital quantities. However, there are three 
difficulties in doing this, all of which relate to the gauge 
choice. First, the T-driver shift condition, which has proven 
very robust and is routinely employed in moving -puncture 
simulations, is usually initialized with a vanishing shift vec- 
tor. Consequently, it takes typically at least one orbit until 
initial transients have decayed in the orbital tracks. It is thus 
not possible, for example, to compute the orbital frequency at 
t = and the required fits can be performed only after the sys- 
tem has reached quasi-equilibrium. Second, depending on the 



choice of a free parameter tj in the T-driver shift, gauge de- 
pendent oscillations at about twice the orbital frequency con- 
taminate the orbital tracks throughout the inspiral for eccen- 
tricities on the order of 10~ 3 and below, which prevents us 
from accurately performing the required fits. Third, while fit- 
ting to time-derivatives simplifies the model functions, it also 
amplifies numerical noise and increases the need for filtering 
even at eccentricities on the order of 5 x 10~ 3 . 

To be specific, we consider one variant of the latest version 
of the method [14], which only makes use of the frequency 
and therefore could in principle be used with a filtered GW 
frequency instead of the orbital frequency to avoid the prob- 
lems mentioned above. This variant requires the following fit 



HnrO ) = S a (t) +B n cos (coat + (j>a) , 



(7.1) 



with fitting parameters Ba, (Oa, and (pa and the following non- 
oscillatory model 
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with k = 1 or k = 2 and free parameters A^ and T c . Various 
choices for the updating formulas for r and Q. at t = can be 
made, but, as far as we are aware, only the following choice 
allows us to rely exclusively on frequency information 



Ar= — — cosfc, 

ZkiQ 

A£l= -— — sinfc- 
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(7.3) 
(7.4) 



When using the time-derivative of the cleaned GW sig- 
nal, one problem with this method is that the GW signal for 
moving -puncture evolutions is extremely noisy for the first 
one to two orbital periods. This implies that the cleaned GW 
signal and the residual eccentricity oscillations will only be 
accurate in the actual fitting window, starting after roughly 
two orbits. The phase offset (j)a of the eccentricity oscilla- 
tions, which is the crucial quantity for this variant of the CCC- 
method, and is determined by the fit at t = 0, may therefore not 
be accurate enough to efficiently reduce the residual eccen- 
tricity. Moreover the frequency (Oa is not constant over time, 
which exacerbates the error in the phase offset further. This 
can be partially addressed by adding a term proportional to t 2 
in the argument of the cosine in ( |7.1| > as mentioned in lfT4l . 

We briefly compare the essence of the variant of the CCC- 
method discussed above to our PN-based eccentricity reduc- 
tion method. The objective of both methods is to compute ap- 
proximate corrections to initial parameters at a starting time 
/ = and a radial separation ro. To do that, residual eccentric- 
ity information from one NR frequency quantity is extracted. 
In the CCC-method this is done by fitting the frequency to a 
nonlinear model function ( |7.1[ ), while in our method the initial 
parameters of a PN model are adjusted such that the residual 
of the time-evolution of this PN-model matches with the NR 
frequency residual. As a second step, updating formulas for 
the initial parameters are used in the CCC-method that rely 
on Newtonian information. In our method, we simply apply 
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the matching perturbation in the opposite direction to recover 
our updated initial parameters. While the computation of the 
match in our method is more involved that the fit in the CCC- 
method, we can avoid a nonlinear fit in the method proper and 
obtain a robust method. 

While there are technical problems in applying the CCC 
method directly to moving-puncture simulations, on the other 
hand it should be possible to apply the method that we have 
presented to SpEC CTS-excision simulations. It would be in- 
teresting in the future to compare the two methods in that set- 
ting. 



VIII. SUMMARY AND DISCUSSION 

We have developed a robust iterative procedure to reduce 
the eccentricity in moving-puncture simulations of black-hole 
binaries. In this method the eccentricity is measured from the 
phase of the GW signal, reducing the gauge dependence of 
previous methods, which used the orbital motion. 

The method relies on calculating differences (residuals) be- 
tween the GW frequency of an NR simulation, and that of an 
analytic model of the same physical system. The key require- 
ments of the model are that it (a) accurately capture the phase 
evolution of the binary; (b) be parametrized by the same initial 
parameters as the NR simulation, which in our case are the ra- 
dial and tangential momenta (p r ,p r ); and (c) we can produce 
a model solution with zero eccentricity. We use solutions of 
a set of PN/EOB equations as the model (Appendix [A}. We 
calculate two residuals with respect to the zero-eccentricity 
model frequency, C0M(p^,Pt)- In one we use the NR GW fre- 
quency, and in the other we use a perturbed model frequency, 
(OM{Kp^,^rP°)- The essence of the method is in finding the 
scale factors (X r ,Xt) such that the two residuals agree in both 
amplitude and phase. The amplitude and phase can be ad- 
justed independently, making it possible to semi-automate this 
procedure. We then update the NR initial parameters with the 
inverse parameters, p r — > p r /X r , p, — > p t /K- We find that with 
this procedure the eccentricity can typically be reduced from 
e i-v 0.01 to e < 10~ 3 in two iteration steps. 

For the method to work, we must filter the NR GW signal, 
and account for dephasing effects in the frequency residuals. 
Even then, noise in the NR waveform prevents us from reduc- 
ing the eccentricity to lower than about 2 x 10~ 4 . However, 
in studies of the mismatch between hybrid PN+NR wave- 
forms, we have seen no evidence that eccentricities as high 
as 0.01 in the final ^10 orbits will have any noticeable in- 
fluence on GW searches or parameter estimation in the Ad- 
vanced detector era. This is somewhat surprising, since at this 
level the eccentricity is visible by eye in the waveform. To 
be conservative we prefer to lower the eccentricity to a level 
where the eccentricity-induced oscillations and secular drifts 
in the GW phase are well below the numerical phase errors in 
our simulations. We choose a tolerance of e ~ 10~ 3 , which 
produces oscillations in the GW phase with an amplitude of 
A<p ~ 0.01 rad during inspiral, and an accumulated phase off- 
set through merger and ringdown of less than 0.2 rad. This is 
well within our numerical phase errors. We note that our anal- 



ysis considered only one configuration, but we do not expect 
the relationship between phase oscillations and eccentricity to 
change very much across the parameter space. 

We have shown that, for typical gauge choices in moving- 
puncture simulations, the frequency of the orbital motion 
cannot be used to accurately measure the eccentricity below 
e ~ 0.002, and even at this level the eccentricity estimator is 
contaminated by gauge effects; see, for example, Fig. [7] 

In large studies of the black-hole-binary parameter space, 
we estimate that the computational overhead in performing 
this eccentricity reduction is between 25% and 35%. 

In our implementation of the method the GW signal is given 
by the Newman-Penrose scalar ^4. We make the important 
observation that the eccentricity measured from the phase of 
W4 will be different to that from the GW strain h. In Ap- 



pendix C 2 we find a simple expression for the relative scale 
factor between the two measures. 

So far the method has been applied to non-precessing bina- 
ries. Precession will introduce additional oscillations into the 
GW amplitude and frequency, which will make it difficult to 
isolate the effects of eccentricity and of precession, and hence 
complicate our method. However, techniques that simplify the 
phase evolution of the precessing-binary waveform, like that 
suggested in [54-56], may alleviate this problem. We will 
consider this further in future work. 

Because our method is applied to the GW signal, it can be 
adapted to any evolution method, and is not limited to moving- 
puncture simulations. It could also be adapted to other com- 
pact binary simulations, for example neutron-star (NS-NS) bi- 
naries, and black-hole-neutron-star (BH-NS) binaries. 
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Appendix A: EOB/PN equations of motion 

For our post-Newtonian evolutions we are using the Hamil- 
tonian equations of motion in the standard Taylor-expanded 
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form, as we have done previously J8] [TOl . and in the EOB 
form [15 16 j. For the Taylor-expanded version, we use the 
non-spinning 3PN accurate Hamiltonian [ 17, 57 58 1 (see also 
I159H6I1 ') and 3.5PN accurate radiation flux ll62"H64l . We add 
both leading-order 1 16 65-69] and next-to-leading order [70- 
l72l contributions to the spin-orbit and spin-spin Hamiltonians, 
and spin-induced radiation flux terms as described in [73 1 (see 
also [68 69 1). In addition we include the flux contribution due 
to the energy flowing in to the black holes, which appears at 
the relative 2.5PN order, as derived in Ref. f74l . 

For the EOB equations, we use the 3PN accurate resum- 
mation [75 1 of the above non-spinning Hamiltonian, and add 
the Taylor-expanded spinning terms. We evolve the resulting 
EOB Hamiltonian evolution equations in the same ADM-TT 
coordinate system used for the Taylor-expanded Hamiltonian 
approach, i.e., representing the EOB momenta and general- 
ized coordinates in terms of the ADM-TT expressions. For 
simplicity, we perform this canonical coordinate transforma- 
tion from the EOB to the ADM-TT phase space only to 2PN 
order, as described in [ 15 1. This is sufficient for our purposes, 
as demonstrated by the successful of our eccentricity reduc- 
tion procedure. In the future we do however plan to use the 
3PN-order transformation [75 1. 



added nonconservative radiation reaction force are 

dH N 
dpi ' 
dH N 



9i = 



dq t 



-F U 



(B2) 
(B3) 



where the "radiation reaction force" is proportional to the mo- 
mentum vector and the rate of energy loss, 



1 dE 
1 QL dt 

The quadrupole flux is given as 

^ = _^ v 2 (Mn) 10/3 

dt 5 



(B4) 



(B5) 



where €l = is the orbital frequency. The radial momentum 
is p r = vMr, and the angular momentum L = p^ is conserved 
for vanishing energy flux. Instead of L we use the tangential 
momentum p, = L/r in the following sections. 



1. Newtonian eccentricity estimators 



Appendix B: Eccentricity perturbations in Newtonian dynamics 

We consider the effect of momenta perturbations on ec- 
centricity in the simple setting of Newtonian binaries. We 
consider both conservative motion, and the inclusion of 
quadrupole radiation reaction. This allows us to derive some 
basic results that were useful in developing and implementing 
our general eccentricity-reduction procedure. 

In Sec. |B l| we first define the customary Newtonian eccen- 
tricity estimators that are based on the analytical solution to 
the Kepler problem linearized in eccentricity. We consider the 
effects of momenta perturbations to conservative dynamics in 
Sec. |B2| and show that the effects of radial and tangential mo- 
menta perturbations are out of phase. We include quadrupole 



radiation reaction in Sec. B 3 which allows us to study the 



dephasing due to eccentricity. 

An understanding of the behavior of GW-signal-based ec- 
centricity estimators requires us to go beyond Newtonian or- 



der, and in Appendix C 1 we consider 1PN effects. This calcu- 
lation allows us to derive the leading-order ratio between the 
eccentricity measured from the GW strain h and the Newman- 
Penrose scalar ^4. 

As a simple model for an inspiraling binary we consider 
Newtonian dynamics with quadrupole radiation reaction as 
discussed in f73l . The Hamiltonian for the Kepler problem 
in polar coordinates q, — (r, 0) in the center-of-mass frame 
with total mass M = m\ +ni2 and symmetric mass-ratio v = 
mxmijM- is 



H\ 



Pr 



Pi 



vM z 



(Bl) 



2vM 2r 2 vM r ' 
where we have set G = 1 . The equations of motion with an 



In this section we collect the customary definitions of ec- 
centricity estimators used in the context of BH binary evolu- 
tions. These estimators are based on Taylor expansions of the 
Kepler solution to linear order in eccentricity. This conserva- 
tive Newtonian setting is the only one in which eccentricity 
can be uniquely defined. 

We first summarize the explicit solutions (see e.g., f76l ) for 
to the Kepler problem. Here we drop the subscript or b for or- 
bital quantities and only give it for the eccentricity estimators. 

The radial separation as a function of the phase is given by 



r(0) 



a(l- 



1 +ecos(p ' 



(B6) 



where a = L 2 /(M 3 v 2 (l — e 2 )) is the semi-major axis, and the 
orbital frequency is 



n(*) 



vMr(<j>) 7 



The Newtonian orbital phase 

(j) =<j) +A e (u), 
is defined in terms of the true anomaly 

/1^^ 1/2 
A e (u) = 2arctan 



1 



tan- 



(B7) 



(B8) 



(B9) 



and eccentric anomaly u, which is related to t by the Kepler 
equation 



j3 := D.Qt = u — esinu, 



(BIO) 



where /3 is the mean anomaly and the average orbital fre- 
quency is Q.q = M/a 3 . 
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Expanding to linear order in the eccentricity, we find 



r(0) 
£2(0) 



L 2 
M 3 v 2 
M 5 V 3 



(l-ecos</>) + ^(e 2 ), 
(l+2ecos0) + ^(e 2 ), 



0(0 =0o + i2 o f + 2esin(QoO + ^(e ), 



and can define eccentricity estimators from the orbital separa- 
tion r, orbital frequency Q. and orbital phase as follows, 
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2. Conservative dynamics 

In this section we study at an analytic level the effect on 
eccentricity of momenta perturbations in a Newtonian binary. 
We start with conservative dynamics. 

For the Kepler problem the orbital eccentricity can be writ- 
ten explicitly as a function of the separation r and the radial 
and tangential momenta p r and p, 



e(j>t,Pr,r) = \ I 



2EL 2 
v 3 M 5 ' 



(B17) 



Note that for a bound system the energy E is negative. 

Circular orbits satisfy p r = and p r = at all times, and 
this requirement leads to the solution p° = vM^M/tq. We 
now perturb either p r or p,, which will give rise to eccen- 
tricity. The eccentricity e(p r ,p t ) has a single minimum at 
the circular-orbit momenta values, which implies that given 
generic elliptical data both momenta need to be adjusted to 
find the circular initial values. 

Fig. [BJillustrates the effect of perturbing the initial tangen- 
tial momentum p t q in a q = 2 binary, where either p t o > p°, 
(left panel), or p t ,o < p° (right panel). Large perturbations 
were chosen to exaggerate the effect. We see that an increase 
in the initial tangential momentum leads to a larger radial sep- 
aration, with a maximum at the phase = n (apastron), while 
the radius keeps the same value at the periastron. In contrast, 
the perturbation of the orbital frequency starts at a maximum 
and reaches a minimum at the apastron. The frequency is also 
shifted by a constant offset. Choosing the initial radial mo- 
mentum smaller than its circular value zero leads to the per- 
turbation of the radial separation that vanish at the periastron 
and apastron and take extremal values half-way in between. 
The perturbation of the orbital frequency behaves similarly, 



but with the opposite sign; see Fig. 15 



For completeness, we also consider the effect of perturb- 
ing the initial separation away from its circular orbit value r°, 
which will also introduce eccentricity. Choosing r^ < r° re- 
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FIG. 14: Newtonian orbits for perturbation of initial tangential mo- 
mentum from circular initial parameters for q = 2,tq = \1M. The 
circular orbit is shown in blue. For the left plot p, q = 1.3/?° and for 
the right plot p, o = 0.7 p° , 
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FIG. 15: Newtonian orbits for perturbation of initial radial momen- 
tum from circular initial parameters for q = 2, r^ = 12M. The circular 
orbit is shown in blue. For the left plot p r Q = — 0.03M and for the 
right plot p r0 = Q.Q3M. 



The perturbation is therefore qualitatively similar to increas- 
ing the tangential momentum from its circular value and the 
associated mode for the orbital frequency is proportional to a 
cosine without the offset that is present in the latter case. This 
is not very surprising, as we still have p r o — and only p t o 
does not have its correct circular value for the new initial ra- 
dial separation. Choosing ro > r° again flips the sign of the 
mode. 
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FIG. 16: Newtonian orbits for a perturbation of the initial radial sep- 
aration from its circular-orbit value for q = 2,r° = 12M. The circular 
orbit is shown in with a solid line. We refer to the radial separation 
of the reference circular parameters as r° . For the left plot r\) — 3r° 
and for the right plot 7'o = 0.6r°. 

We can analyze the analytical solution for the Kepler prob- 



sults in an elliptic orbit that starts at the periastron (see Fig 16 1, lem given in section B 1 to gain more insight into the behavior 
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of these momenta perturbations. For our purposes it is suffi- 
cient to do this up to linear order in eccentricity. To simplify 
the solutions we consider the circular initial momenta for a 
given initial separation ro 



l.Orr 



p°=0, 

p° = VMy/Mfo, 



(B18) 
(B19) 



as the reference values and write the reference frequency as 
£1° := QlPt>P° = M 5 v 3 /L° 3 . Perturbing the initial tangential 
momentum p t p at t = by a factor A, is equivalent to perturb- 
ing the angular momentum L by the same factor, as long as the 
initial separation stays fixed. If we choose X, > 1 and 0o = 0, 
the orbital frequency perturbation is 



s?^' 1 ■=q a -' p ' -a 



pAtpf ,_c\%t 

Q. 
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2e(p°,X,p°,r ) 



COS 0(f) 



A t / Af 

D.° (A, - 1) (-3 +4cos0(f)) + 0{{Xt ~ I) 2 )- 



(B20) 



Similarly, we can perturb the radial momentum p r from its 
circular value zero. Since p° = 0, we cannot scale it to obtain 
an eccentric solution. Therefore, we choose p r Q < and find 



= 2Q.°e(p r fl,p°,r ) sin <j)(t) 



Z ) 



(B21) 



:2Q.°\p r \ 



'o 



v 2 M 3 



sin<H0 + ^ 2 )- 



Note that Eqn. ( |B6[ ) assumes that the motion starts at the peri- 
astron. Choosing p r a < and 0o = the periastron is shifted 
by 7T/2, hence the sine. 

The perturbation of p, is associated with a cosine mode 
with an offset, while the perturbation of p r gives rise to a sine 



mode. In Fig. 17 we plot frequency residuals for perturbations 
of p t , p r , and r, which show this mode behavior. This is con- 
sistent with the heuristic argument presented at the beginning 



of Sec. IV A This fact allows us to extract the oscillatory in- 



formation of a single frequency quantity to adjust two physical 
parameters (i.e., the radial and tangential momenta), and is a 
key feature of our eccentricity reduction method. One could 
also develop a simpler method, based on the variation of only 
one parameter (for example the initial binary separation), but 
such a method would not be able to reduce the eccentricity 
between some bound. See Sec. [VIA] where we estimate that 
this bound is around e ~ 0.005. 



3. Newtonian evolutions with quadrupole flux 

In this section we include quadrupole radiation-reaction. 
This leads to dephasing between the quasi-circular and eccen- 
tric configurations. 

The addition of a radiation reaction through Eqns. |B4} and 
to the Newtonian equations of motion (|B2|i leads to in- 
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FIG. 17: Frequency residuals for perturbations of the tangential mo- 
mentum, radial momentum, and radial separation for a conservative 
Newtonian model. The residuals are shown for a.q = 2,ro= 12M bi- 
nary and perturbations p t q = 1 .03/^ (red curve), p r Q = —0.01 (green 
long-dashed curve) and rg = 0.9r° (blue short-dashed curve). 



spiral, and the radial separation for a quasi-circular (QC) con- 
figurations is 



'vM 3 

,- ( >, = 4(^L(r c -f; 



1/4 



(B22) 



(see e.g., 11771 1. where T c is the coalescence time and can be 
expressed in terms of the initial radial separation as 



^ = 250^3-' 



(B23) 



We choose the initial tangential momentum as in the conser- 
vative case, by Eqn. ( |B19[ ). The initial radial momentum that 
will result in QC inspiral is now nonzero and can be found 
by combining equations ( |B22| i and ( |B23[ ) and taking a time 
derivative 



Pr.QC = vMr (t = Q;T = T c (r )) = - 



64M' 
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In Fig. [18] we show the radial separation and the orbital 
frequency for quasi-circular (QC) inspiral and perturbations 
thereof, which lead to eccentric inspiral. We show the fre- 
quency residuals for these cases in Fig. 19 (top panel). It is 



apparent that it is only the perturbation of the tangential mo- 
mentum that causes a significant 'dephasing' (red line). For 
the simple Newtonian-plus-quadrupole-flux model, a pertur- 
bation of p r causes no dephasing in the frequency and a per- 
turbation of r only a very slight dephasing. For perturbations 
of PN/EOB approxmiants, however, significant dephasing is 
generic and therefore we need to correct for it. The simple 
model considered here serves mainly as an illustration of the 
dephasing effect. 

We can calculate the dephasing effect in our Newtonian 
model in the following way. First we find the orbital frequency 
for QC inspiral by combining Kepler's law Q? = M/P with 
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FIG. 18: Radial separation (top) and orbital frequency (bottom) for 
Newtonian plus quadrupole flux inspiral for QC data (q = 2,rr, = 
12M) and perturbations p,Q = 1.01p,,Qc, Pr.o = 5p,-,QC an d nj = 
0.98rQc from the QC values. 



the evolution of the radial separation ( |B22[ > to obtain 
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Writing the coalescence time as in ( |B23| >, it is apparent that the 
single initial parameter ro in the orbital frequency ( |B25[ ) can 
be perturbed. (Alternatively, T c could we written as a func- 
tion of initial frequency, so that CIq can be perturbed.) We ex- 
pand the dephasing 2>a and @n = Q.(t; T c (r 8)) - Cl(f, T c (r )) 
around 8 = 1 and obtain 
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1 3/8 



whereC=i5 3 / 8 [l/(vM 5 /3)] ; 

The functional form of the dephasing can be used as a fitting 
model beyond its Newtonian setting by letting the prefactors 
of the rational terms be free parameters. We then arrive at the 



following model 



%2,model(f;<5,r c ,A,B) 



+ 



A(8-l)T c 

"(r c -r)"/8 

B{8-l) 2 (6tT c + 5T c 2 ) 



(T c -t 



19/8 



(B27) 



with parameters 8, T C ,A,B. 

We can use this rational model to remove the dephasing, as 
shown in the bottom panel of Fig. 19 However, when apply- 



ing our eccentricity reduction method to full NR simulations, 
we find that nonlinear fits to this model are difficult to tune, 
and in practice we have found that the polynomial fit ( |2.19[ ) 
to capture the dephasing behavior just as well, while also be- 
ing more robust. The only drawback of the polynomial fit is 
that it results in artifacts at the beginning and end of the fitting 
window, where the polynomial is prone to picking up parts of 
oscillations in the residuals. This problem can be alleviated 
by discarding part of the fitting window near the boundary. 
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FIG. 19: Frequency residuals for perturbations of tangential mo- 
mentum, radial momentum, and radial separation for the Newtonian- 
plus-quadrupole-flux model. The residuals are computed from per- 
turbations p, o = l.Olp, qc, Pr.O = 5Pr,QC ani ^ r = 0.98r<jc from 
QC data for a q = 2, ro = V2M binary. The upper panel illustrates the 
dephasing effect. In the lower panel the dephasing has been removed 
by a fit to lfB27]>. 
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Appendix C: Eccentricity for 1PN compact binaries 

This discussion is based on the 1PN equations of motion 
in a quasi-Keplerian parametrization as given in [78 1, with 
the notation in same places modified for consistency with the 
rest of this paper. After defining the relevant quantities we 
linearize in the "time-eccentricity" <?, and evaluate the usual 
Newtonian eccentricity estimators with the orbital phase and 
frequency. Instead of e t we could have chosen the "radial ec- 
centricity" e r or the "phase eccentricity" e§ defined below. 

In section |C2| we use the quadrupole formula to obtain the 
strain polarizations h + and h x from the 1PN accurate orbital 
phase. Furthermore, we calculate the phases of h and ¥4 to 
linear order in <?, , and evaluate the usual eccentricity estimator 
for the GW phase. Finally, we relate the results of the orbital 
and GW eccentricity estimators for the 1PN compact binary as 
functions of the periastron advance parameter k. This result is 
surprising at first glance: if we measure the eccentricity from 
the GW phase, we obtain different answers if we use the strain 
h and the Newman-Penrose scalar v ?4. It is important to bear 
this result in mind whenever the GW signal is used to measure 
the eccentricity (as for example in [7|); fortunately to leading 
order the two can be related. 



1. Definitions 

The following is based on the treatment in IT781 . 
We define the orbital phase 

<l> rb = <h + (l+k)A e , (CI) 

the true anomaly A e in terms of the eccentric anomaly u 

'11 \ 1/2 
1 +e± \ ' 



A e (u) = 2arctan 



1-, 



tan- 



and the mean anomaly 

j3 :=n(t — to) = u — e t smu. 
The radial separation is given by 

r = a r (1 — e,-cos«) , 
where 



3/ kv-1 

a r = H 

k\ 6 4 



(C2) 



(C3) 



(C4) 



(C5) 



The eccentricities e&, e r and e t can be related in terms of the 
periastron advance parameter k 

^ = l + jkl-«?)(-2(v-4)), (C6) 

e t o 

- = 1 + ^(1 - e , 2 )(8-3v). (C7) 

e t o 

The frequency of radial eccentricity oscillations is 

n = a r = 2n/P r , (C8) 



where P r measures the time between two consecutive perias- 
tron passages, while the average orbital frequency is 



fy = (l+it)Q r . 



(C9) 



The fractional periastron advance per orbit, k, depends on the 
frequency of eccentricity oscillations 



k = 3e z 



» 2 / 3 



(CIO) 



where e counts the order of the inverse speed of light c, e = 
l/cl 

Expanding up to linear order in eccentricity (the different 
eccentricities can be related to each other) we find for the or- 
bital phase 

<j> orb = fo + iy + e,(l +k) (l + ^) sin(£V). (Cll) 



The Newtonian definition of the orbital phase eccentricity 
then yields 

Kb- Corbie = 0) 1+e^/e, 
e<j>,orb= z =e t (\+k) . (C12) 

The orbital frequency is simply the time-derivative of the or- 
bital phase 



n = %+e,^ 1 + — cos(O-0 (C13) 



and the associated eccentricity is 

Q-Q(e = 0) l + e^/e, 



en = 



2Q 



= e,- 



(C14) 



The ratio between the orbital phase and frequency eccentrici- 
ties is 



^(j),orb 

en 



= l+k, 



(C15) 



as given in [7|. At Newtonian order there is no periastron 
advance, and k = 0, and, as we have already emphasized, all 
eccentricity estimators agree. 



2. Eccentricity estimators for the GW signal 

In this section we calculate the ratio between the eccen- 
tricities measured from the phase of the GW strain h and the 
Newman-Penrose scalar v ?4. 

We first calculate expressions for ¥4 and the strain h 



by combining the 1PN phase given in Sec. C 1 with the 
quadrupole formula 



,rr = ^r ( ,_ r)j 



(C16) 



where Z, ; is the reduced quadrupole moment and TT projects 
out the transverse traceless part of a tensor. The strain can be 
decomposed as 



hff =h + e±+h x efj, 



(C17) 
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in terms of two polarization tensors ef- and e* 11771 . 

For an observer in the wave zone at a distance r and incli- 
nation angle 9 relative to the plane in which the binary orbits, 
the two polarization modes of the strain are 
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sin(2 f)(-M + r 2 <p-r 2> 



-2rr0cos(20) \ (C19) 

The only spherical-harmonic modes of interest are the 
(£,m) — (2, ±2) modes. The (£,m) — (2,±1) modes van- 
ish irrespective of the choice of 9 at 1PN order, and the 
(£,m) = (2,0) mode is real and does not contribute to the 
phase. To make the computation of the phases tractable we as- 
sume that the binary is optimally oriented to the observer and 
set 0=0, which implies that only the (£,m) = (2,2) mode 
remains. This is sufficient, because out eccentricity reduction 
method considers only the frequency from ^4,22 ■ 

The wave strain 

h = h+-ih x (C20) 

is related to the Newman-Penrose scalar W4 by 

¥ 4 = h + -ih x . (C21) 

The phase of ^4 and h can be obtained locally in time by 
computing the complex argument of either W4 or h. 

Applying the customary definition of an eccentricity esti- 
mator for the GW phase 



e<i>,GW 



0GW — <Pfit 



we find 



and 



2</>[fc] = e f 



3+4/fce 2 



3 m~. 



21 + lOke 
' 16" 



2 



(C22) 



(C23) 



(C24) 



We may Taylor expand the factor between e^iyj and e^i 
in £ up to lPN-order (e 2 = 1/c 2 ) to obtain an approximation 
in the low k or low frequency limit 



? 0[V4] 

e<j,[h] 



1 
4 



3ke 2 



(C25) 
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FIG. 20: Comparison of eccentricity estimators e*ri»>i and e*iu for 



the reference (q = 1,%\ = 0,^2 = 0.25) configuration with e ■■ 
The theoretical scale factor between them is 1.36 (see text). 



0.01. 



It is clear that in general we must specify whether the GW- 
phase eccentricity was measured from h or ^4. In the main 
text we refer only to the 'GW phase', 0gw> by which we mean 
the phase computed from W4. 

Figure 20 compares the eccentricity estimators e^yA an d 
e^r/,1 for our reference example q = 2,%\ = 0,^2 = configu- 
ration, with high eccentricity e = 0.01. We take k — 0.4 from 
the NR results in Fig. 7 of [7 1 (assuming it still holds for q = 2) 
for an average orbital frequency in the early evolution of 
Mco ~ 0.027. We then calculate the ratio e^^j/e^/,] ~ 1.36. 
This factor agrees with the numerical data to about 15% accu- 
racy in the amplitude. We have repeated the analysis for the 
same configuration with eccentricity e = 0.006, and find sim- 
ilar agreement. For smaller eccentricities the noise in the ^4 
residual makes the comparison less conclusive. 

For the sake of completeness, we compute the eccentricity 
estimator for the frequency of v ?4, 



?ea[ift] 



4< 21 



like 1 



(C26) 



Up to 1PN order the ratio between esuu ,] and e a \«, ■.] equals 
1 + k, as expected. 

In Sec. [V] we compare the eccentricity measured from both 
1 , 4 and the orbital-motion frequency Q.. The ratio between 
the phase- and frequency-based eccentricities is l+k, as dis- 
cussed earlier in Appendix C 1 But we are now considering 



the phase and frequency from different physical quantities, re- 
spectively ^4 and the orbital motion, and so we do not expect 
that the same relationship will hold. 

To obtain an approximation for the low-£ or low-frequency 
limit, we Taylor expand the factor betwee n e^ W i \ and <?q up to 
lPN-order in e. We combine Eqns. (|C6]), ( [CT4| , and ( |Cl4| , 



k:= 



2 t>{¥4. 

en 



21 
16 



7V - 1 ]^ 

32 4, 



(C27) 



The average orbital frequency for the NR q = 1 example is 
Mco fa 0.027. We can recover k at the 1PN level by combining 



23 



Eqns. (|C8|)-(|C10|l to give 



k = 3 



0.027 \ 2/3 



(C28) 



which yields k w 0.23. Then, K « 1.27. Taking fc from the 
Fig. 7 of Ref. Q, a value k ~ 0.38 seems reasonable. This 



gives a factor K k, 1 .25, a bit smaller than the factor 1 + k. 

We find good agreement with a factor 1 .25 between eccen- 
tricities in Figs|9]and 10 The disagreement is at most 10%, 
which is inside the error bars of the eccentricity estimators 
discussed in section [VTB Moreover, we expect that the 1PN 
expressions also contribute a significant error. 
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